home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 26 / Cream of the Crop 26.iso / os2 / octa209s.zip / octave-2.09 / src / xpow.cc < prev    next >
C/C++ Source or Header  |  1996-10-12  |  15KB  |  767 lines
  1. /*
  2.  
  3. Copyright (C) 1996 John W. Eaton
  4.  
  5. This file is part of Octave.
  6.  
  7. Octave is free software; you can redistribute it and/or modify it
  8. under the terms of the GNU General Public License as published by the
  9. Free Software Foundation; either version 2, or (at your option) any
  10. later version.
  11.  
  12. Octave is distributed in the hope that it will be useful, but WITHOUT
  13. ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  14. FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  15. for more details.
  16.  
  17. You should have received a copy of the GNU General Public License
  18. along with Octave; see the file COPYING.  If not, write to the Free
  19. Software Foundation, 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
  20.  
  21. */
  22.  
  23. #ifdef HAVE_CONFIG_H
  24. #include <config.h>
  25. #endif
  26.  
  27. #include <cassert>
  28. #include <climits>
  29.  
  30. #include "CColVector.h"
  31. #include "CDiagMatrix.h"
  32. #include "CMatrix.h"
  33. #include "EIG.h"
  34. #include "dDiagMatrix.h"
  35. #include "dMatrix.h"
  36. #include "oct-cmplx.h"
  37.  
  38. #include "error.h"
  39. #include "ov.h"
  40. #include "utils.h"
  41. #include "xpow.h"
  42.  
  43. static inline int
  44. xisint (double x)
  45. {
  46.   return (D_NINT (x) == x
  47.       && ((x >= 0 && x < INT_MAX)
  48.           || (x <= 0 && x > INT_MIN)));
  49. }
  50.  
  51. // Safer pow functions.
  52. //
  53. //       op2 \ op1:   s   m   cs   cm
  54. //            +--   +---+---+----+----+
  55. //   scalar   |     | 1 | 5 |  7 | 11 |
  56. //                  +---+---+----+----+
  57. //   matrix         | 2 | * |  8 |  * |
  58. //                  +---+---+----+----+
  59. //   complex_scalar | 3 | 6 |  9 | 12 |
  60. //                  +---+---+----+----+
  61. //   complex_matrix | 4 | * | 10 |  * |
  62. //                  +---+---+----+----+
  63.  
  64. // -*- 1 -*-
  65. octave_value
  66. xpow (double a, double b)
  67. {
  68.   if (a < 0.0 && (int) b != b)
  69.     {
  70.       Complex atmp (a);
  71.       return pow (atmp, b);
  72.     }
  73.   else
  74.     return pow (a, b);
  75. }
  76.  
  77. // -*- 2 -*-
  78. octave_value
  79. xpow (double a, const Matrix& b)
  80. {
  81.   octave_value retval;
  82.  
  83.   int nr = b.rows ();
  84.   int nc = b.cols ();
  85.  
  86.   if (nr == 0 || nc == 0 || nr != nc)
  87.     error ("for x^A, A must be square");
  88.   else
  89.     {
  90.       EIG b_eig (b);
  91.       ComplexColumnVector lambda (b_eig.eigenvalues ());
  92.       ComplexMatrix Q (b_eig.eigenvectors ());
  93.  
  94.       for (int i = 0; i < nr; i++)
  95.     {
  96.       Complex elt = lambda (i);
  97.       if (imag (elt) == 0.0)
  98.         lambda (i) = pow (a, real (elt));
  99.       else
  100.         lambda (i) = pow (a, elt);
  101.     }
  102.       ComplexDiagMatrix D (lambda);
  103.  
  104.       retval = ComplexMatrix (Q * D * Q.inverse ());
  105.     }
  106.  
  107.   return retval;
  108. }
  109.  
  110. // -*- 3 -*-
  111. octave_value
  112. xpow (double a, const Complex& b)
  113. {
  114.   Complex result;
  115.   Complex atmp (a);
  116.   result = pow (atmp, b);
  117.   return result;
  118. }
  119.  
  120. // -*- 4 -*-
  121. octave_value
  122. xpow (double a, const ComplexMatrix& b)
  123. {
  124.   octave_value retval;
  125.  
  126.   int nr = b.rows ();
  127.   int nc = b.cols ();
  128.  
  129.   if (nr == 0 || nc == 0 || nr != nc)
  130.     error ("for x^A, A must be square");
  131.   else
  132.     {
  133.       EIG b_eig (b);
  134.       ComplexColumnVector lambda (b_eig.eigenvalues ());
  135.       ComplexMatrix Q (b_eig.eigenvectors ());
  136.  
  137.       for (int i = 0; i < nr; i++)
  138.     {
  139.       Complex elt = lambda (i);
  140.       if (imag (elt) == 0.0)
  141.         lambda (i) = pow (a, real (elt));
  142.       else
  143.         lambda (i) = pow (a, elt);
  144.     }
  145.       ComplexDiagMatrix D (lambda);
  146.  
  147.       retval = ComplexMatrix (Q * D * Q.inverse ());
  148.     }
  149.  
  150.   return retval;
  151. }
  152.  
  153. // -*- 5 -*-
  154. octave_value
  155. xpow (const Matrix& a, double b)
  156. {
  157.   octave_value retval;
  158.  
  159.   int nr = a.rows ();
  160.   int nc = a.cols ();
  161.  
  162.   if (nr == 0 || nc == 0 || nr != nc)
  163.     error ("for A^b, A must be square");
  164.   else
  165.     {
  166.       if ((int) b == b)
  167.     {
  168.       int btmp = (int) b;
  169.       if (btmp == 0)
  170.         {
  171.           retval = DiagMatrix (nr, nr, 1.0);
  172.         }
  173.       else
  174.         {
  175.           // Too much copying?
  176.           // XXX FIXME XXX -- we shouldn't do this if the exponent is
  177.           // large...
  178.  
  179.           Matrix atmp;
  180.           if (btmp < 0)
  181.         {
  182.           btmp = -btmp;
  183.  
  184.           int info;
  185.           double rcond = 0.0;
  186.  
  187.           atmp = a.inverse (info, rcond, 1);
  188.  
  189.           if (info == -1)
  190.             warning ("inverse: matrix singular to machine\
  191.  precision, rcond = %g", rcond);
  192.         }
  193.           else
  194.         atmp = a;
  195.  
  196.           Matrix result (atmp);
  197.           for (int i = 1; i < btmp; i++)
  198.         result = result * atmp;
  199.  
  200.           retval = result;
  201.         }
  202.     }
  203.       else
  204.     {
  205.       EIG a_eig (a);
  206.       ComplexColumnVector lambda (a_eig.eigenvalues ());
  207.       ComplexMatrix Q (a_eig.eigenvectors ());
  208.  
  209.       for (int i = 0; i < nr; i++)
  210.         lambda (i) = pow (lambda (i), b);
  211.  
  212.       ComplexDiagMatrix D (lambda);
  213.  
  214.       retval = ComplexMatrix (Q * D * Q.inverse ());
  215.     }
  216.     }
  217.  
  218.   return retval;
  219. }
  220.  
  221. // -*- 6 -*-
  222. octave_value
  223. xpow (const Matrix& a, const Complex& b)
  224. {
  225.   octave_value retval;
  226.  
  227.   int nr = a.rows ();
  228.   int nc = a.cols ();
  229.  
  230.   if (nr == 0 || nc == 0 || nr != nc)
  231.     error ("for A^b, A must be square");
  232.   else
  233.     {
  234.       EIG a_eig (a);
  235.       ComplexColumnVector lambda (a_eig.eigenvalues ());
  236.       ComplexMatrix Q (a_eig.eigenvectors ());
  237.  
  238.       for (int i = 0; i < nr; i++)
  239.     lambda (i) = pow (lambda (i), b);
  240.  
  241.       ComplexDiagMatrix D (lambda);
  242.  
  243.       retval = ComplexMatrix (Q * D * Q.inverse ());
  244.     }
  245.  
  246.   return retval;
  247. }
  248.  
  249. // -*- 7 -*-
  250. octave_value
  251. xpow (const Complex& a, double b)
  252. {
  253.   Complex result;
  254.  
  255.   if (xisint (b))
  256.     result = pow (a, (int) b);
  257.   else
  258.     result = pow (a, b);
  259.  
  260.   return result;
  261. }
  262.  
  263. // -*- 8 -*-
  264. octave_value
  265. xpow (const Complex& a, const Matrix& b)
  266. {
  267.   octave_value retval;
  268.  
  269.   int nr = b.rows ();
  270.   int nc = b.cols ();
  271.  
  272.   if (nr == 0 || nc == 0 || nr != nc)
  273.     error ("for x^A, A must be square");
  274.   else
  275.     {
  276.       EIG b_eig (b);
  277.       ComplexColumnVector lambda (b_eig.eigenvalues ());
  278.       ComplexMatrix Q (b_eig.eigenvectors ());
  279.  
  280.       for (int i = 0; i < nr; i++)
  281.     {
  282.       Complex elt = lambda (i);
  283.       if (imag (elt) == 0.0)
  284.         lambda (i) = pow (a, real (elt));
  285.       else
  286.         lambda (i) = pow (a, elt);
  287.     }
  288.       ComplexDiagMatrix D (lambda);
  289.  
  290.       retval = ComplexMatrix (Q * D * Q.inverse ());
  291.     }
  292.  
  293.   return retval;
  294. }
  295.  
  296. // -*- 9 -*-
  297. octave_value
  298. xpow (const Complex& a, const Complex& b)
  299. {
  300.   Complex result;
  301.   result = pow (a, b);
  302.   return result;
  303. }
  304.  
  305. // -*- 10 -*-
  306. octave_value
  307. xpow (const Complex& a, const ComplexMatrix& b)
  308. {
  309.   octave_value retval;
  310.  
  311.   int nr = b.rows ();
  312.   int nc = b.cols ();
  313.  
  314.   if (nr == 0 || nc == 0 || nr != nc)
  315.     error ("for x^A, A must be square");
  316.   else
  317.     {
  318.       EIG b_eig (b);
  319.       ComplexColumnVector lambda (b_eig.eigenvalues ());
  320.       ComplexMatrix Q (b_eig.eigenvectors ());
  321.  
  322.       for (int i = 0; i < nr; i++)
  323.     {
  324.       Complex elt = lambda (i);
  325.       if (imag (elt) == 0.0)
  326.         lambda (i) = pow (a, real (elt));
  327.       else
  328.         lambda (i) = pow (a, elt);
  329.     }
  330.       ComplexDiagMatrix D (lambda);
  331.  
  332.       retval = ComplexMatrix (Q * D * Q.inverse ());
  333.     }
  334.  
  335.   return retval;
  336. }
  337.  
  338. // -*- 11 -*-
  339. octave_value
  340. xpow (const ComplexMatrix& a, double b)
  341. {
  342.   octave_value retval;
  343.  
  344.   int nr = a.rows ();
  345.   int nc = a.cols ();
  346.  
  347.   if (nr == 0 || nc == 0 || nr != nc)
  348.     error ("for A^b, A must be square");
  349.   else
  350.     {
  351.       if ((int) b == b)
  352.     {
  353.       int btmp = (int) b;
  354.       if (btmp == 0)
  355.         {
  356.           retval = DiagMatrix (nr, nr, 1.0);
  357.         }
  358.       else
  359.         {
  360.           // Too much copying?
  361.           // XXX FIXME XXX -- we shouldn't do this if the exponent is
  362.           // large...
  363.  
  364.           ComplexMatrix atmp;
  365.           if (btmp < 0)
  366.         {
  367.           btmp = -btmp;
  368.  
  369.           int info;
  370.           double rcond = 0.0;
  371.  
  372.           atmp = a.inverse (info, rcond, 1);
  373.  
  374.           if (info == -1)
  375.             warning ("inverse: matrix singular to machine\
  376.  precision, rcond = %g", rcond);
  377.         }
  378.           else
  379.         atmp = a;
  380.  
  381.           ComplexMatrix result (atmp);
  382.           for (int i = 1; i < btmp; i++)
  383.         result = result * atmp;
  384.  
  385.           retval = result;
  386.         }
  387.     }
  388.       else
  389.     {
  390.       EIG a_eig (a);
  391.       ComplexColumnVector lambda (a_eig.eigenvalues ());
  392.       ComplexMatrix Q (a_eig.eigenvectors ());
  393.  
  394.       for (int i = 0; i < nr; i++)
  395.         lambda (i) = pow (lambda (i), b);
  396.  
  397.       ComplexDiagMatrix D (lambda);
  398.  
  399.       retval = ComplexMatrix (Q * D * Q.inverse ());
  400.     }
  401.     }
  402.  
  403.   return retval;
  404. }
  405.  
  406. // -*- 12 -*-
  407. octave_value
  408. xpow (const ComplexMatrix& a, const Complex& b)
  409. {
  410.   octave_value retval;
  411.  
  412.   int nr = a.rows ();
  413.   int nc = a.cols ();
  414.  
  415.   if (nr == 0 || nc == 0 || nr != nc)
  416.     error ("for A^b, A must be square");
  417.   else
  418.     {
  419.       EIG a_eig (a);
  420.       ComplexColumnVector lambda (a_eig.eigenvalues ());
  421.       ComplexMatrix Q (a_eig.eigenvectors ());
  422.  
  423.       for (int i = 0; i < nr; i++)
  424.     lambda (i) = pow (lambda (i), b);
  425.  
  426.       ComplexDiagMatrix D (lambda);
  427.  
  428.       retval = ComplexMatrix (Q * D * Q.inverse ());
  429.     }
  430.  
  431.   return retval;
  432. }
  433.  
  434. // Safer pow functions that work elementwise for matrices.
  435. //
  436. //       op2 \ op1:   s   m   cs   cm
  437. //            +--   +---+---+----+----+
  438. //   scalar   |     | * | 3 |  * |  9 |
  439. //                  +---+---+----+----+
  440. //   matrix         | 1 | 4 |  7 | 10 |
  441. //                  +---+---+----+----+
  442. //   complex_scalar | * | 5 |  * | 11 |
  443. //                  +---+---+----+----+
  444. //   complex_matrix | 2 | 6 |  8 | 12 |
  445. //                  +---+---+----+----+
  446. //
  447. //   * -> not needed.
  448.  
  449. // -*- 1 -*-
  450. octave_value
  451. elem_xpow (double a, const Matrix& b)
  452. {
  453.   octave_value retval;
  454.  
  455.   int nr = b.rows ();
  456.   int nc = b.cols ();
  457.  
  458.   // For now, assume the worst.
  459.  
  460.   if (a < 0.0)
  461.     {
  462.       Complex atmp (a);
  463.       ComplexMatrix result (nr, nc);
  464.       for (int j = 0; j < nc; j++)
  465.     for (int i = 0; i < nr; i++)
  466.       result (i, j) = pow (atmp, b (i, j));
  467.  
  468.       retval = result;
  469.     }
  470.   else
  471.     {
  472.       Matrix result (nr, nc);
  473.       for (int j = 0; j < nc; j++)
  474.     for (int i = 0; i < nr; i++)
  475.       result (i, j) = pow (a, b (i, j)); 
  476.  
  477.       retval = result;
  478.     }
  479.  
  480.   return retval;
  481. }
  482.  
  483. // -*- 2 -*-
  484. octave_value
  485. elem_xpow (double a, const ComplexMatrix& b)
  486. {
  487.   int nr = b.rows ();
  488.   int nc = b.cols ();
  489.  
  490.   ComplexMatrix result (nr, nc);
  491.   for (int j = 0; j < nc; j++)
  492.     for (int i = 0; i < nr; i++)
  493.       result (i, j) = pow (a, b (i, j));
  494.  
  495.   return result;
  496. }
  497.  
  498. // -*- 3 -*-
  499. octave_value
  500. elem_xpow (const Matrix& a, double b)
  501. {
  502.   octave_value retval;
  503.  
  504.   int nr = a.rows ();
  505.   int nc = a.cols ();
  506.  
  507.   if ((int) b != b && a.any_element_is_negative ())
  508.     {
  509.       ComplexMatrix result (nr, nc);
  510.       for (int j = 0; j < nc; j++)
  511.     for (int i = 0; i < nr; i++)
  512.       {
  513.         Complex atmp (a (i, j));
  514.         result (i, j) = pow (atmp, b);
  515.       }
  516.  
  517.       retval = result;
  518.     }
  519.   else
  520.     {
  521.       Matrix result (nr, nc);
  522.       for (int j = 0; j < nc; j++)
  523.     for (int i = 0; i < nr; i++)
  524.       result (i, j) = pow (a (i, j), b);
  525.  
  526.       retval = result;
  527.     }
  528.  
  529.   return retval;
  530. }
  531.  
  532. // -*- 4 -*-
  533. octave_value
  534. elem_xpow (const Matrix& a, const Matrix& b)
  535. {
  536.   octave_value retval;
  537.  
  538.   int nr = a.rows ();
  539.   int nc = a.cols ();
  540.  
  541.   int b_nr = b.rows ();
  542.   int b_nc = b.cols ();
  543.  
  544.   if (nr != b_nr || nc != b_nc)
  545.     {
  546.       gripe_nonconformant ("operator .^", nr, nc, b_nr, b_nc);
  547.       return octave_value ();
  548.     }
  549.  
  550.   int convert_to_complex = 0;
  551.   for (int j = 0; j < nc; j++)
  552.     for (int i = 0; i < nr; i++)
  553.       {
  554.     double atmp = a (i, j);
  555.     double btmp = b (i, j);
  556.     if (atmp < 0.0 && (int) btmp != btmp)
  557.       {
  558.         convert_to_complex = 1;
  559.         goto done;
  560.       }
  561.       }
  562.  
  563. done:
  564.  
  565.   if (convert_to_complex)
  566.     {
  567.       ComplexMatrix complex_result (nr, nc);
  568.  
  569.       for (int j = 0; j < nc; j++)
  570.     for (int i = 0; i < nr; i++)
  571.       {
  572.         Complex atmp (a (i, j));
  573.         Complex btmp (b (i, j));
  574.         complex_result (i, j) = pow (atmp, btmp);
  575.       }
  576.  
  577.       retval = complex_result;
  578.     }
  579.   else
  580.     {
  581.       Matrix result (nr, nc);
  582.  
  583.       for (int j = 0; j < nc; j++)
  584.     for (int i = 0; i < nr; i++)
  585.       result (i, j) = pow (a (i, j), b (i, j));
  586.  
  587.       retval = result;
  588.     }
  589.  
  590.   return retval;
  591. }
  592.  
  593. // -*- 5 -*-
  594. octave_value
  595. elem_xpow (const Matrix& a, const Complex& b)
  596. {
  597.   int nr = a.rows ();
  598.   int nc = a.cols ();
  599.  
  600.   ComplexMatrix result (nr, nc);
  601.   for (int j = 0; j < nc; j++)
  602.     for (int i = 0; i < nr; i++)
  603.       result (i, j) = pow (a (i, j), b);
  604.  
  605.   return result;
  606. }
  607.  
  608. // -*- 6 -*-
  609. octave_value
  610. elem_xpow (const Matrix& a, const ComplexMatrix& b)
  611. {
  612.   int nr = a.rows ();
  613.   int nc = a.cols ();
  614.  
  615.   int b_nr = b.rows ();
  616.   int b_nc = b.cols ();
  617.  
  618.   if (nr != b_nr || nc != b_nc)
  619.     {
  620.       gripe_nonconformant ("operator .^", nr, nc, b_nr, b_nc);
  621.       return octave_value ();
  622.     }
  623.  
  624.   ComplexMatrix result (nr, nc);
  625.   for (int j = 0; j < nc; j++)
  626.     for (int i = 0; i < nr; i++)
  627.       result (i, j) = pow (a (i, j), b (i, j));
  628.  
  629.   return result;
  630. }
  631.  
  632. // -*- 7 -*-
  633. octave_value
  634. elem_xpow (const Complex& a, const Matrix& b)
  635. {
  636.   int nr = b.rows ();
  637.   int nc = b.cols ();
  638.  
  639.   ComplexMatrix result (nr, nc);
  640.   for (int j = 0; j < nc; j++)
  641.     for (int i = 0; i < nr; i++)
  642.       {
  643.     double btmp = b (i, j);
  644.     if (xisint (btmp))
  645.       result (i, j) = pow (a, (int) btmp);
  646.     else
  647.       result (i, j) = pow (a, btmp);
  648.       }
  649.  
  650.   return result;
  651. }
  652.  
  653. // -*- 8 -*-
  654. octave_value
  655. elem_xpow (const Complex& a, const ComplexMatrix& b)
  656. {
  657.   int nr = b.rows ();
  658.   int nc = b.cols ();
  659.  
  660.   ComplexMatrix result (nr, nc);
  661.   for (int j = 0; j < nc; j++)
  662.     for (int i = 0; i < nr; i++)
  663.       result (i, j) = pow (a, b (i, j));
  664.  
  665.   return result;
  666. }
  667.  
  668. // -*- 9 -*-
  669. octave_value
  670. elem_xpow (const ComplexMatrix& a, double b)
  671. {
  672.   int nr = a.rows ();
  673.   int nc = a.cols ();
  674.  
  675.   ComplexMatrix result (nr, nc);
  676.  
  677.   if (xisint (b))
  678.     {
  679.       for (int j = 0; j < nc; j++)
  680.     for (int i = 0; i < nr; i++)
  681.       result (i, j) = pow (a (i, j), (int) b);
  682.     }
  683.   else
  684.     {
  685.       for (int j = 0; j < nc; j++)
  686.     for (int i = 0; i < nr; i++)
  687.       result (i, j) = pow (a (i, j), b);
  688.     }
  689.  
  690.   return result;
  691. }
  692.  
  693. // -*- 10 -*-
  694. octave_value
  695. elem_xpow (const ComplexMatrix& a, const Matrix& b)
  696. {
  697.   int nr = a.rows ();
  698.   int nc = a.cols ();
  699.  
  700.   int b_nr = b.rows ();
  701.   int b_nc = b.cols ();
  702.  
  703.   if (nr != b_nr || nc != b_nc)
  704.     {
  705.       gripe_nonconformant ("operator .^", nr, nc, b_nr, b_nc);
  706.       return octave_value ();
  707.     }
  708.  
  709.   ComplexMatrix result (nr, nc);
  710.   for (int j = 0; j < nc; j++)
  711.     for (int i = 0; i < nr; i++)
  712.       {
  713.     double btmp = b (i, j);
  714.     if (xisint (btmp))
  715.       result (i, j) = pow (a (i, j), (int) btmp);
  716.     else
  717.       result (i, j) = pow (a (i, j), btmp);
  718.       }
  719.  
  720.   return result;
  721. }
  722.  
  723. // -*- 11 -*-
  724. octave_value
  725. elem_xpow (const ComplexMatrix& a, const Complex& b)
  726. {
  727.   int nr = a.rows ();
  728.   int nc = a.cols ();
  729.  
  730.   ComplexMatrix result (nr, nc);
  731.   for (int j = 0; j < nc; j++)
  732.     for (int i = 0; i < nr; i++)
  733.       result (i, j) = pow (a (i, j), b);
  734.  
  735.   return result;
  736. }
  737.  
  738. // -*- 12 -*-
  739. octave_value
  740. elem_xpow (const ComplexMatrix& a, const ComplexMatrix& b)
  741. {
  742.   int nr = a.rows ();
  743.   int nc = a.cols ();
  744.  
  745.   int b_nr = b.rows ();
  746.   int b_nc = b.cols ();
  747.  
  748.   if (nr != b_nr || nc != b_nc)
  749.     {
  750.       gripe_nonconformant ("operator .^", nr, nc, b_nr, b_nc);
  751.       return octave_value ();
  752.     }
  753.  
  754.   ComplexMatrix result (nr, nc);
  755.   for (int j = 0; j < nc; j++)
  756.     for (int i = 0; i < nr; i++)
  757.       result (i, j) = pow (a (i, j), b (i, j));
  758.  
  759.   return result;
  760. }
  761.  
  762. /*
  763. ;;; Local Variables: ***
  764. ;;; mode: C++ ***
  765. ;;; End: ***
  766. */
  767.